somehow the tests are failing tho... I'll figure it out

Yeah, it looks like that function is actually being called with some kind of array of vectors, so it needs to be multidimension-aware?

OK, I see the issue...

It's that quadilateral_area gives angle_between 3D arrays... (e.g. arrays of np.shape = (3, nx, ny)) while angle_between expects vectors of size 3.

How can I make it broadcast over the 2nd and 3rd dimension? I don't remember how to do it in python... :(

I guess we could do np.cross(..., axis=0) and np.tensordot(..., axis=([0],[0])), but that's a bit ugly... It might also work if we use arrays of shape (nx, ny, 3) instead? That would also be the "correct" indexing convention for numpy -- having the vector axis first is the Fortran convention.

This

    nv, nx, ny = np.shape(c0)

    c0 = np.reshape(c0, (nx, ny, nv))
    c1 = np.reshape(c1, (nx, ny, nv))
    c2 = np.reshape(c2, (nx, ny, nv))
    c3 = np.reshape(c3, (nx, ny, nv))

reshapes the inputs... But then how do I ask the angle_between to broadcast over last two dimensions?

I think if you write it in such a way that the dot products and cross products are computed over the last axis, it'll just work as expected. I think that's the case for the cross product, from https://numpy.org/doc/stable/reference/generated/numpy.cross.html#numpy.cross:

The cross product of a and b in is a vector perpendicular to both a and b. If a and b are arrays of vectors, the vectors are defined by the last axis of a and b by default, and these axes can have dimensions 2 or 3.

But I'm not sure what to do for the dot product...

So at the moment quadrilateral_area computes an array of areas.

I suggest we change it so that

def quadilateral_area(v1, v2, v3, v4):
    """Returns area of a spherical quadrilateral on the unit sphere that
    has vertices on 3-vectors `v1`, `v2`, `v3`, `v4` (counter-clockwise
    orientation is implied). The area is computed as the excess of the
    sum of the spherical angles of the quadrilateral from 2π."""

    a1 = angle_between(v1, v2, v4)
    a2 = angle_between(v2, v3, v1)
    a3 = angle_between(v3, v4, v2)
    a4 = angle_between(v4, v1, v3)

    return a1 + a2 + a3 + a4 - 2.0 * np.pi

that gives the area of a single quadrilateral. Then we have

def quadilateral_areas(lat, lon):
    ....

that uses quadilateral_area to computes the quadrilateral areas for the grid.

How do you like that?

OK, I did it but it's not using broadcasting but it's doing the dreaded for-loop. See

@angus-g or @ashjbarnes can you convert it to broadcasting or something more efficient?

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.39% 🎉

Comparison is base (2f9d12b) 15.42% compared to head (8d45e19) 15.81%.

@@            Coverage Diff             @@
##             main      #65      +/-   ##
==========================================
+ Coverage   15.42%   15.81%   +0.39%     
==========================================
  Files           2        3       +1     
  Lines         428      430       +2     
==========================================
+ Hits           66       68       +2     
  Misses        362      362              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Yeah, I'll take a look.

@angus-g I added a validation in quadilateral_area via d091f32 that checks if all vectors provided are of the same length. Can we add a test to ensure that en exception is raised when vectors of different length are provided?

@angus-g I added a validation in quadilateral_area via d091f32 that checks if all vectors provided are of the same length. Can we add a test to ensure that en exception is raised when vectors of different length are provided?

chatGPT told me how to do it! done via aa5ed7a

I like it - very verbose and readable

OK, just the optimization of the for-loop with broadcasting and we are ready to merge then :)

Here's an attempt at doing the broadcasting: dccdede

It's a little ugly to have to dstack everything in the tests, but on a 50x50 grid it takes 1.6ms/it vs. 715ms/it (400x speedup) vs. the for-loop.